//  -Wl,--gc-sections
//  -ffunction-sections
#include "ncbi/serial/iterator.hpp"
#include "ncbi/objmgr/object_manager.hpp"
#include "ncbi/objmgr/scope.hpp"
#include "ncbi/objmgr/util/sequence.hpp"
#include "ncbi/objmgr/seq_vector.hpp"
#include "ncbi/algo/blast/api/bl2seq.hpp"
#include "ncbi/algo/blast/api/blast_options.hpp"
#include "ncbi/objects/seq/Bioseq.hpp"
#include "ncbi/objects/seq/NCBIeaa.hpp"
#include "ncbi/objects/seq/Seq_inst.hpp"
#include "ncbi/objects/seq/Seqdesc.hpp"
#include "ncbi/objects/seq/Seq_descr.hpp"
#include "ncbi/objects/seqalign/Seq_align.hpp"
#include "ncbi/objects/seqalign/Dense_seg.hpp"
#include "ncbi/objects/seqalign/Score.hpp"
#include "ncbi/objects/general/Object_id.hpp"
#include "ncbi/util/sequtil/sequtil_convert.hpp"
#include "ncbi/objects/seq/IUPACna.hpp"
#include "ncbi/objects/seq/NCBI4na.hpp"
#include "ncbi/objects/seq/NCBI2na.hpp"


#include "bl2seq_mine.h"


#include <ctype.h>
#include <boost/foreach.hpp>
#include <string>
#include <algorithm>
#include <stdexcept>
#include <iostream>
#include <vector>
# define DEBUG_LEVEL 1
#include <debug_func.h>
using std::cerr;
using std::string;
using std::vector;
using std::max;
using namespace ncbi::blast;
using namespace ncbi::objects;
using namespace ncbi;
using bcpp_align::t_align_results;
using bcpp_align::t_gap;

// requires the following NCBI libraries
//-lxblast -lxnetblast -lseqdb -lblastdb -lscoremat -lxobjutil -ltables -lxconnect -lxcompress -lxobjmgr -lseqset -lseq -lseqcode -lsequtil -lpub -lmedline -lbiblio -lgeneral -lxser -lxutil -lxncbi -ldl


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
// t_blast_options
//
//

// pimpl to avoid any NCBI dependencies in main code
struct t_blast_options::t_impl
{
	public:
		static CRef<CObjectManager>     m_ObjMgr;  ///< instance of the object manager
		CRef<CBlastOptionsHandle>       m_options_handle;
};


CRef<CObjectManager> t_blast_options::t_impl::m_ObjMgr;

t_blast_options::t_blast_options(   bool seg_filter ,
									double e_value_cutoff ,
									const string& matrix,
									double theoretical_db_size,
									unsigned gap_x_dropoff,
									unsigned wordsize,
									unsigned gap_open,
									unsigned gap_extend)
	: p_impl(new t_impl())
{
	debug_func_cerr( 8 );

	p_impl->m_ObjMgr = CObjectManager::GetInstance();
	if (!p_impl->m_ObjMgr)
		throw std::runtime_error("Could not initialize blast object manager");

	if (matrix != "DNA")
	{
		p_impl->m_options_handle.Reset(CBlastOptionsFactory::Create(ncbi::blast::eBlastp));
		DNA = false;
	}
	else
	{
		debug_cerr( 8, "Dna mode" );
		p_impl->m_options_handle.Reset(CBlastOptionsFactory::Create(ncbi::blast::eBlastn));
		DNA = true;
	}
	if ( !p_impl->m_options_handle )
		throw std::runtime_error("Could not create blast options");

	CBlastOptions& opt = p_impl->m_options_handle->SetOptions();
	if (!is_dna())
	{
		debug_cerr( 8, matrix.c_str() << " mode" );
		opt.SetMatrixName(matrix.c_str());
	}
	else
	{
		gap_open = 2;
		gap_extend = 2;
		//opt.SetMismatchPenalty(-1);
		opt.SetMatchReward(2);
	}

	opt.SetFilterString(seg_filter ? "T" : "F");
	debug_cerr( 8, "seg_filter? " << (seg_filter ? "T" : "F"));
	debug_cerr( 8, matrix);

	// optional
	if (gap_open)
	{
		opt.SetGapOpeningCost(gap_open);
	}
	if (gap_extend)
	{
		opt.SetGapExtensionCost(gap_extend);
	}
	if (wordsize)
	{
		opt.SetWordSize(wordsize);
	}

	//opt.SetGapExtnAlgorithm(eGreedyWithTracebackExt);
	//opt.SetUngappedExtension(false);


	opt.SetGapXDropoffFinal(gap_x_dropoff);
	opt.SetGapXDropoff(gap_x_dropoff);
	opt.SetXDropoff(gap_x_dropoff);
	debug_cerr( 8, "gap_x_dropoff=" << gap_x_dropoff);
	opt.SetEvalueThreshold(e_value_cutoff);
	debug_cerr( 8, "e_value_cutoff=" << e_value_cutoff);

	//opt.SetEffectiveSearchSpace(theoretical_db_size);
	opt.SetDbLength(theoretical_db_size);
	debug_cerr( 8, "theoretical_db_size=" << theoretical_db_size);



	opt.SetGappedMode(true);


	/* Immediate greedy gapped extension with traceback */
		//opt.SetGapExtnAlgorithm(eGreedyWithTracebackExt);
//		opt.SetUngappedExtension(false);
	/* Two-step greedy extension, no ungapped extension */
//		opt.SetGapExtnAlgorithm(eGreedyExt);
//		opt.SetGapTracebackAlgorithm(eGreedyTbck);
//		opt.SetUngappedExtension(false);
//	/* Two-step greedy extension after ungapped extension*/
//		opt.SetGapExtnAlgorithm(eGreedyExt);
//		opt.SetGapTracebackAlgorithm(eGreedyTbck);

};



//
// t_blast_options
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888














//________________________________________________________________________________________


//	get_seqlocvec

//		create a NCBI TSeqLocVector using sequence / accession

//________________________________________________________________________________________
namespace 
{
	
	TSeqLocVector
	get_seqlocvec(CObjectManager& objmgr, bool dna, const string& accession, const string& sequence)
	{
		debug_func_cerr( 8 );

			
		CRef<CSeq_entry> seq_entry(new CSeq_entry);
		CRef<CBioseq>          seq(new CBioseq);
		{
			seq->SetInst().SetRepr(CSeq_inst::eRepr_raw);
			if (!dna)
			{
				seq->SetInst().SetMol(CSeq_inst::eMol_aa);
				seq->SetInst().SetSeq_data().SetNcbieaa().Set(sequence);
				seq->SetInst().SetLength(sequence.length());
			}
			else
			{
				seq->SetInst().SetMol(CSeq_inst::eMol_dna);
				//seq->SetInst().SetSeq_data().SetIupacna().Set(sequence);

				vector<char> v((sequence.length() + 1) / 2, '\0');
				CSeqUtil::ECoding coding;
				CSeqConvert::Pack(sequence, 
								  CSeqUtil::e_Iupacna, v, coding);
				if (coding == CSeqUtil::e_Ncbi2na) 
				{
					seq->SetInst().SetSeq_data().SetNcbi2na().Set().assign(	v.begin(),
																			v.begin() + (sequence.length() + 3) / 4);
//					seq->SetInst().SetLength((sequence.length() + 3) / 4);
					seq->SetInst().SetLength(sequence.length());
				}
				else
				{
					swap(seq->SetInst().SetSeq_data().SetNcbi4na().Set(), v);
					seq->SetInst().SetLength(sequence.length());
				}
				//seq->SetInst().SetSeq_data().SetNcbieaa().Set(sequence);
			}
			seq->SetId().push_back(CRef<CSeq_id>
						  (new CSeq_id(CSeq_id::e_Local, accession, kEmptyStr)));
	
			seq_entry->SetSeq(*seq);
		}
	
		CRef<CScope> scope(new CScope(objmgr));
		scope->AddDefaults();
		scope->AddTopLevelSeqEntry(*seq_entry);

		CRef<CSeq_loc> seqloc(new CSeq_loc());
		seqloc->SetInt().SetTo(sequence.length() - 1);
		seqloc->SetInt().SetFrom(0);
		seqloc->SetInt().SetId().Assign(*seq_entry->GetSeq().GetId().front());
		if (dna)
		{
			seqloc->SetInt().SetStrand(eNa_strand_plus);
		}
		else
		{		
			seqloc->SetInt().SetStrand(eNa_strand_unknown); //eNa_strand_unknown
		}
		TSeqLocVector retval;
		retval.push_back(SSeqLoc(seqloc, scope));
	
		return retval;
	}
	

	bool retrieve_data(TSeqAlignVector 		seqalignv,
						t_align_results&	results,
						unsigned 			seq_len1,
						unsigned 			seq_len2)
	{
		debug_func_cerr( 8 );
		if (!seqalignv.size())
			throw std::runtime_error("Bl2Seq internal error: No SeqAlignVector");
		const CSeq_align_set& seq_aligns = *seqalignv[0];
		if (!seq_aligns.Get().size())
			return false;
		CRef<CSeq_align> hit_align = seq_aligns.Get().front();
		const CSeq_align_set& disc = hit_align->GetSegs().GetDisc();
		if (!disc.Get().size())
			return false;
		CRef<CSeq_align> seqalign = disc.Get().front();
	
		// get dense segment representing top hit
		const CDense_seg& dense_seg = seqalign->GetSegs().GetDenseg();
	
		// get scores
		unsigned cnt_identical = std::numeric_limits<unsigned>::max();
		results.e__value = -9999;
		results.bitscore = -9999;
		BOOST_FOREACH(CRef<CScore> score, seqalign->GetScore())
		{
			if (score->GetId().GetStr() == "e_value")
				results.e__value = score->GetValue().GetReal();
			else if (score->GetId().GetStr() == "bit_score")
				results.bitscore = score->GetValue().GetReal();
			else if (score->GetId().GetStr() == "num_ident")
				cnt_identical = score->GetValue().GetInt();
		}
		assert(cnt_identical != std::numeric_limits<unsigned>::max());
		assert(results.e__value != -9999);
		assert(results.bitscore != -9999);
	
		// split starts into two
		const std::vector< TSignedSeqPos >& starts = dense_seg.GetStarts();
		std::vector< TSignedSeqPos > starts1;
		std::vector< TSignedSeqPos > starts2;
		debug_cerr( 8, "Starts");
		for (unsigned i = 0; i < starts.size(); i += 2)
		{
			debug_cerr( 8, i << "," << starts[i] << "," << starts[i + 1] << "\n");
			starts1.push_back(starts[i]);
			starts2.push_back(starts[i + 1]);
		}
			
		const std::vector< TSeqPos >& lens = dense_seg.GetLens();
		// no alignment
		if (!lens.size())
			throw std::runtime_error("Bl2Seq internal error: No alignment segment lengths");
		for (unsigned i = 0; i < lens.size(); i++)
			debug_cerr( 8, i << "," << lens[i] << "\n");
		
	
		results.align_len = 0;
		unsigned ungapped_len1 = seq_len1;
		unsigned ungapped_len2 = seq_len2;
	
	
		// length of current gaps
		unsigned gap_len1 = 0;
		unsigned gap_len2 = 0;
		results.align_beg = max(starts1[0], starts2[0]);
		if (starts2[0] > starts1[0])
			gap_len1 = starts2[0] - starts1[0];
		else
			gap_len2 = starts1[0] - starts2[0];
	
		// go through segments, counting gaps
		for (unsigned i = 0; i < lens.size(); i++)
		{
	
			results.align_len += lens[i];
	
			// save gaps
			if (starts1[i] == -1)
				gap_len1 = lens[i];
			else if (starts2[i] == -1)
				gap_len2 = lens[i];
	
			// start of non-gap position: append saved gap parameter to string
			else
			{
				results.aligned_residues += lens[i];
				if (gap_len1)
				{
					results.gaps1.push_back(t_gap(starts1[i], gap_len1));
					debug_cerr(10, "Gaps1 Start = " << starts1[i] << ", len = " << gap_len1);
					seq_len1 += gap_len1;
					debug_cerr(10, "seq len1 = " << seq_len1);
				}
				else if (gap_len2)
				{
					results.gaps2.push_back(t_gap(starts2[i], gap_len2));
					debug_cerr(10, "Gaps2 Start = " << starts2[i] << ", len = " << gap_len2);
					seq_len2 += gap_len2;
					debug_cerr(10, "seq len2 = " << seq_len2);
				}
				gap_len1 = gap_len2 = 0;
			}
		}
	
		// get beginning gap
		unsigned beg_gap_len = 0;
		if (results.gaps1.size() && results.gaps1[0].beg == 0)
			beg_gap_len = results.gaps1[0].len;
		else if (results.gaps2.size() && results.gaps2[0].beg == 0)
			beg_gap_len = results.gaps2[0].len;
		debug_cerr(10, "beg_gap_len = " << beg_gap_len);
	
		results.total_len = results.align_beg + results.align_len + abs((int)seq_len2 - (int)seq_len1);
		debug_cerr(10, "total_len = " << results.total_len);
		debug_cerr(10, "align_len = " << results.align_len);
		debug_cerr(10, "seq_len1 = " << seq_len1);
		debug_cerr(10, "seq_len2 = " << seq_len2);
	
		// add final gap
		unsigned end_gap_len = 0;
		if (seq_len1 < seq_len2)
		{
			end_gap_len = seq_len2 - seq_len1;
			results.gaps1.push_back(t_gap(ungapped_len1, end_gap_len));
			debug_cerr(10, "Gaps1 Start = " << ungapped_len1 << ", len = " << end_gap_len << "\n");
		}
		else if (seq_len1 > seq_len2)
		{
			end_gap_len = seq_len1 - seq_len2;
			results.gaps2.push_back(t_gap(ungapped_len2, end_gap_len));
			debug_cerr(10, "Gaps2 Start = " << ungapped_len2 << ", len = " << end_gap_len << "\n");
		}
		debug_cerr(10, "end_gap_len = " << end_gap_len);
		results.len1 = ungapped_len1;
		results.len2 = ungapped_len2;
		debug_cerr(10, "ungapped_len1 = " << ungapped_len1);
		debug_cerr(10, "ungapped_len2 = " << ungapped_len2);
	
		results.percent_id = static_cast<double>(cnt_identical) / results.align_len;
		return true;
	}
	
	
	
}










//________________________________________________________________________________________


//	do_bl2seq


//________________________________________________________________________________________
bool do_bl2seq (t_blast_options& blast_options, const string& seq1, const string& seq2,
				t_align_results& results)
{
	debug_func_cerr( 8 );
	// Retrieve input sequences
	TSeqLocVector query_loc = get_seqlocvec(*blast_options.p_impl->m_ObjMgr, 
											blast_options.is_dna(),  "A", seq1);
	TSeqLocVector subject_loc = get_seqlocvec(*blast_options.p_impl->m_ObjMgr, 
												blast_options.is_dna(),  "B", seq2);
	//string n_seq_1 = seq1;
	//string n_seq_2 = seq2;
	//bcpp_align::bl_prepare_sequence_str(n_seq_1);
	//bcpp_align::bl_prepare_sequence_str(n_seq_2);
	//TSeqLocVector query_loc = get_seqlocvec(*blast_options.p_impl->m_ObjMgr, 
	//										blast_options.is_dna(),  "A", n_seq_1);
	//TSeqLocVector subject_loc = get_seqlocvec(*blast_options.p_impl->m_ObjMgr, 
	//											blast_options.is_dna(),  "B", n_seq_2);
	debug_cerr( 8, "Create blaster..." );
	CBl2Seq blaster(query_loc, subject_loc, blast_options.p_impl->m_options_handle.GetObject());
	TSeqAlignVector seqalignv;
	try
	{
		debug_cerr( 8, "Run blaster..." );
		seqalignv = blaster.Run();
		debug_cerr( 8, "Blaster finished..." );
	}
	catch(std::exception& e)
	{
		// strip error of NCBI preamble and rethrow
		string message = e.what();
		string::size_type pos = message.find("Error: ");
		if (pos != string::npos)
			message.erase(0,  pos + 7);
		throw std::runtime_error("?" + message);
	}
	// hack to prevent bl2seq from delete option
	return retrieve_data(seqalignv, results, seq1.length(), seq2.length());
}














/*
TSeqAlignVector
	.push_back(CRef<CSeq_align_set> seq_aligns
		->Set().push_back(
			CRef<CSeq_align> hit_align;
				->SetType(CSeq_align::eType_disc);
				->SetDim(2);         // BLAST only creates pairwise alignments
				->SetSegs().SetDisc().Set().push_back(
					CRef<CSeq_align> seqalign;
						->SetType(CSeq_align::eType_partial);
						->SetDim(2);         // BLAST only creates pairwise alignments
						->SetSegs().SetDenseg(CRef<CDense_seg> dense_seg
							->SetDim(2);                       
							->SetLens() = vector<TSeqPos>& lengths
							->SetStarts() = vector<TSignedSeqPos>& starts
							->SetNumseg((int) lengths.size()); 
						->SetScore().push_back(
							CRef<CScore> score
								->SetId(CRef<CObject_id> id
									->SetStr("e_value/bit_score/score/num_ident");
								->SetValue(CRef<CScore::C_Value> val
									->SetReal(d);


	// Our poor man's formatting ...  
	for (unsigned int index = 0; index < seqalignv.size(); ++index)
	{
		const CSeq_align_set& seq_aligns = *seqalignv[index];
		BOOST_FOREACH(CRef<CSeq_align> hit_align, seq_aligns.Get()) 
		{
			const CSeq_align_set& disc = hit_align->GetSegs().GetDisc();
			BOOST_FOREACH(CRef<CSeq_align> seqalign, disc.Get()) 
			{
				const CDense_seg& dense_seg = seqalign->GetSegs().GetDenseg();
				const std::vector< TSeqPos >& lens = dense_seg.GetLens();
				std::cout << "Lengths=";
				BOOST_FOREACH(const TSeqPos& pos, lens)
					std::cout  <<  pos << ", " ;
				std::cout  << std::endl;
				const std::vector< TSignedSeqPos >& starts = dense_seg.GetStarts();
				std::cout << "Starts=";
				BOOST_FOREACH(const TSignedSeqPos& s, starts)
					std::cout  <<  s << ", " ;
				std::cout  << std::endl;

				BOOST_FOREACH(CRef<CScore> score, seqalign->GetScore())
				{
					std::cout  << score->GetId().GetStr() << " = " ;
					if (score->GetValue().IsReal())
						std::cout << score->GetValue().GetReal()  << std::endl;
					else
						std::cout << score->GetValue().GetInt()  << std::endl;
				}
			}

		}
	}
*/



